【麻将】计算普通立直麻将规则下天和的概率

在普通立直麻将规则下，庄家天和的概率是多少？

随便上网一查资料就知道，这是 $\frac{2143179701279}{708384171528036000}$，约等于 $\frac{1}{330530}$。

那么怎么计算这个东西呢

假如我们枚举庄家的 $14$ 张配牌，然后逐一检测是否和牌并计算概率。不考虑赤宝牌，那么不同的配牌有方案 $326520504500$ 种。这是 $34$ 个 $[0,4]$ 内的整数之和等于 $14$ 的方案数，枚举大于 $4$ 的数的个数容斥，可以得到计算式 $\sum_{t=0}^{2}(-1)^tC(34,t)C(34+14-1-5t,34-1)=326520504500$。这样枚举量过于巨大，在普通的计算机上运行起码需要数小时的时间。

换一种思路，枚举和牌的牌型。国士无双型和牌型有 $13$ 种，七对子型有 $C(34,7)=5379616$ 种，一般型的计算复杂一些，我们可以大致估计一下：$34$ 种雀头，$21$ 种顺子，$34$ 种刻子，我们可以依次枚举雀头和四组面子。面子之间不分顺序，可以规定一个顺序并从小到大枚举，过程中筛去不合法的方案和重复的方案，这样枚举量大约在 $34\times\frac{(34+21)^4}{4!}\approx 1.3\times 10^7$ 这个数量级，实现得当的程序运行数秒钟即可得出结果。实际上一般型有 $11498658$ 种和牌型，其中 $4668$ 种同时也符合七对子型，避免重复计算需要减去。

这 $326520504500$ 种牌型里，每种牌型并不是等概率的。概率可以这样计算：假设所有牌都不一样，那配牌有 $C(136,14)=4250305029168216000$ 种。如果一种牌型里一种牌出现了 $k$ 张，那它的方案数就要乘以 $C(4,k)$。例如，每一种七对子型出现的概率都是 $\frac{C(4,2)^7}{C(136,14)}=\frac{279936}{4250305029168216000}$。

剩下的就是写程序了，其实还是没那么容易写的，代码如下：

#include <cstdio>
#include <unordered_set>
#include <algorithm>
typedef long long ll;
const int c4[5] = {1,4,6,4,1};
int p[4]; ll g[4];
// p=牌型数 c=概率乘以C(136,14)
// 0=国士型 1=七对子型 2=一般型 3=一般型内的七对子型（负值）
int a[34]; // 0~26三色数牌，27~33字牌
ll hsh = 0, h[34];
std::unordered_set<ll> T; // 哈希判重
void find(int d, int x) {
	if (!x) {
		if (T.count(hsh)) return;
		T.insert(hsh);
		ll ng = 1;
		bool ct = true;
		for (int i = 0; i < 34; i++) {
			ct &= (a[i] == 0 || a[i] == 2);
			ng *= c4[a[i]];
		}
		p[2]++;
		g[2] += ng;
		if (ct) {
			p[3]--;
			g[3] -= ng;
		}
		return;
	}
	// 34~54三色顺子
	for (int i = d; i < 55; i++) {
		if (i < 34) {
			if (a[i] <= 1) {
				a[i] += 3;
				hsh += 3 * h[i];
				find(i+1, x-1);
				hsh -= 3 * h[i];
				a[i] -= 3;
			}
		} else {
			int y = (i-34) + (i-34)/7 * 2;
			if (a[y] <= 3 && a[y+1] <= 3 && a[y+2] <= 3) {
				a[y]++; a[y+1]++; a[y+2]++;
				hsh += 31 * h[y];
				find(i, x-1);
				hsh -= 31 * h[y];
				a[y]--; a[y+1]--; a[y+2]--;
			}
		}
	}
}
int main()  {
	p[0] = 13;
	g[0] = p[0] * 100663296ll; // C(4,1)^12 * C(4,2)
	p[1] = 5379616; // C(34,7)
	g[1] = p[1] * 279936ll; // C(4,2)^7
	h[0] = 1;
	for (int i = 1; i < 34; i++)
		h[i] = h[i-1] * 5;
	for (int i = 0; i < 34; i++) {
		a[i] += 2;
		hsh += 2 * h[i];
		find(0, 4);
		hsh -= 2 * h[i];
		a[i] -= 2;
	}
	int ap = 0; ll ag = 0;
	for (int i = 0; i < 4; i++) {
		printf("%d %lld\n", p[i], g[i]);
		ap += p[i]; ag += g[i];
	}
	printf("%d %lld\n", ap, ag);
	ll q = 4250305029168216000; // C(136,14);
	ll gcd = std::__gcd(ag, q);
	ag /= gcd; q /= gcd;
	printf("%lld/%lld\n", ag, q);
}

程序输出：

13 1308622848
5379616 1505948184576
11498658 11353128141498
-4668 -1306741248
16873619 12859078207674
2143179701279/708384171528036000

可以看到输出和上面给出的数据都是一致的，那个 $\frac{2143179701279}{708384171528036000}$ 就是最终结果了